%% Function file that calculates concentration (using backscatter--TSS calibration (regression from TSS calibration), cov stress, TKE stress,
%and settling velocity (filtering out the washload concentration) from the profiler ADV 

%Created by Cartwright, updated by Fall November 2016

function  getAverageADV_cov_using_backConc(filename,backConc,b)

%gets the burst average values for sample heights and calculates the
%bulk settling velocity for the burst

eval (['load ' filename '.mat'])
eval (['load ' filename 'index.asc'])

%b=[-0.877 0.0337];    % For Sontek curve ADV backscatter to TSS calibration using all 2007 + 2008 survey data - set to 0 if no calibration available
%b=[-2.273 0.0447];    % ADV backscatter using Nortek curve (KAF March 2014)

mncorr=70;     % minimum acceptable correlation values - used to eliminate individual burst values that don't meet criteria - set to 0 to do no QC at this step

i = eval ([ filename 'index']);

 sizei=size(i)
 
 count=0;
 
 for ni=3:2:sizei
     count=count+1
     clear data good advdatg gd u v w vel v2 ur2 vr wr ur3 vr2 wr2 t trend
     data(:,1)= Vx(i(ni):i(ni+1)); 
     data(:,2)= Vy(i(ni):i(ni+1)); 
     data(:,3)= Vz(i(ni):i(ni+1));
     data(:,4)= A1(i(ni):i(ni+1));
     data(:,5)= A2(i(ni):i(ni+1)); 
     data(:,6)= A3(i(ni):i(ni+1)); 
     data(:,7)= C1(i(ni):i(ni+1)) ;
     data(:,8)= C2(i(ni):i(ni+1));
     data(:,9)= C3(i(ni):i(ni+1)) ;
     data(:,10)= heading(i(ni):i(ni+1)) ;
     data(:,11)= pitch(i(ni):i(ni+1)); 
     data(:,12)= roll(i(ni):i(ni+1)); 
     data(:,13)= ADVtemp(i(ni):i(ni+1)) ;
     data(:,14)= ADVpress(i(ni):i(ni+1));
    
     [p,q]=size(data);
% use specified minimum correlation to eliminate "bad" data
good= data(:,7)>=mncorr & data(:,8)>=mncorr & data(:,9)>=mncorr;  % indices corresponding to all beams correlation greater than minimum
advdatg=data(good,:);

gd=find(good==1);
nsmp=length(gd);       % number of samples remaining after quality control and used for averaging

if nsmp > 3 
    
% correct velocities for pitch, roll, and heading
%prad=advdatg(:,11).*2.*pi./360;  %pitch in radians
%rrad=advdatg(:,12).*2.*pi./360;  %roll in radians


u=advdatg(:,1);
v=advdatg(:,2);
w=advdatg(:,3);

%new matrix approach developed by Jarrell Smith 10/10/2012

vmean=[mean(u);mean(v);mean(w)]; %[u;v;w] %mean values
[theta,phi,vmag]=cart2sph(vmean(1),vmean(2),vmean(3)); %axis angles
ct=cos(-theta); st=sin(-theta);
cp=cos(phi); sp=sin(phi);
Rz=[ct,-st,0;st,ct,0;0,0,1]; %rotation about z
Ry=[cp,0,sp;0,1,0;-sp,0,cp]; %rotation about y
vmean2=Ry*(Rz*vmean); %test transformation (not required)
%to get timeseries
vel=[u(:)';v(:)';w(:)']; %raw timeseries
v2=Ry*(Rz*vel); %transformed timeseries

ur2=v2(1,:)';
vr=v2(2,:)';
wr=v2(3,:)';

figure
plot(ur2)
hold on
plot(vr,'r-')
plot(wr,'g-')
legend('ur2','vr','wr')
hold off
pause
close all

%this was used by pat dickhudt to rotate velocity data
%[u,v,w]=pitchcorr(u,v,w,prad,rrad);

%[u,v,w]=rollcorr(u,v,w,prad,rrad);

%[thba,sba]=cart2pol(mean(u),mean(v));

%ur=u.*cos(-thba)-v.*sin(-thba);    % rotated u velocity
%vr=u.*sin(-thba)+v.*cos(-thba);    % rotated v velocity

%uba2=mean(ur);
%wba=mean(w);
%[thba2,sba2]=cart2pol(uba2,wba);
%ur2=ur.*cos(-thba2)-w.*sin(-thba2);      % vertically rotated u velocity guarantees u is parallel to bed
%wr=ur.*sin(-thba2)+w.*cos(-thba2);       % vertically rotated w velocity guarantees u is perpendicular to bed

% correct rotated velocities for drift within the burst
t=1:length(ur2);
a=polyfit(t',ur2,1);
trend=a(1)*t' + a(2);
ur3=ur2-trend;

mean_velocity = mean(ur2);  %%% print out

a=polyfit(t',vr,1);
trend=a(1)*t' + a(2);
vr2=vr-trend;

a=polyfit(t',wr,1);
trend=a(1)*t' + a(2);
wr2=wr-trend;

%convert the direction theta (polar) into degree North coordinate.
thetaN =(theta*360./2./pi);
if(thetaN<0);
thetaN = 360 + thetaN;
end

if(thetaN>=0&thetaN<90);
compass_head = abs(thetaN - 90);
end

if(thetaN>=90&thetaN<=360);
compass_head = abs(450 - thetaN);
end

%Combine getBackscatterAvgSD_SE into this function
ampmn=mean(advdatg(:,4:6),2);      % determine mean backscatter amplitude from 3 beams
ampstd_rcd=mean(std(advdatg(:,4:6),1));  %mean of the standard deviation for the 3 beams

bscatter_stats.ampmn=ampmn;                         %burst mean backscatter amplitude
bscatter_stats.std_ampmn=std(ampmn);                %burst standard deviation
bscatter_stats.stderr_ampmn=std(ampmn)/sqrt(nsmp);  %burst standard error
bscatter_stats.ampstd_rcd=ampstd_rcd;               %mean record standard deviation
bscatter_stats.nsmp_good= nsmp;                      %number of good samples
bscatter_stats.nsmp= p ;                             %number of samples in the burst

Cb=exp(b(2).*ampmn+b(1));         % apply calibration to get burst TSS

mean_conc = mean(Cb);  %%%% print out

% subsititute nan matrix for covariance matrix if 1 or less bursts left
% after QC
if nsmp<=3
cv=nan(4,4);
else
cv=cov([Cb ur3 vr2 wr2]);  
end


%Calculate Settling Velocities and save in dateout
backConc2=min(Cb);                         %%% minimum concentration inferred from ADV backscatter
Wspump = 10*cv(1,4)./(mean(Cb)-backConc);  %%%  Ws  directly observed settling velocity using backCONC from TSS pump samples (mm/s)
Wsobs = 10*cv(1,4)./(mean(Cb)-backConc2);  %%% Ws  directly observed settling velocity; back conc is the min CB observed over each depth per burst
Ws_ALL= 10*cv(1,4)./(mean(Cb)); %get settling velocity with background concentration included. Newly added 8/11/14.

% Calculate Reynolds bed stresse
P=1010/100/100; % conversion from u'w' in (cm/s)^2 to Pa (kg/ms2)--need to covernt cm2 to m2 because density in kg/m3
Re_stress=-cv(2,4); %  u'w' directly observed turbulent Reynolds stress in (cm/s)^2
%tau_cov=P*cv(2,4);   % Define directly observed stress in Pascals COV method
tau_cov=P*-Re_stress;   % Think above line should be -cv(2,4) kaf May2014



%% Save data in matrix "datout"
datout(count,1)=mean(ur2);     % burst averaged U velocity
datout(count,2)=mean(vr);     % burst averaged V velocity
datout(count,3)=mean(wr);     % burst averaged W velocity
datout(count,4)=mean(advdatg(:,4));     % burst averaged x beam backscatter amplitude
datout(count,5)=mean(advdatg(:,5));     % burst averaged y beam backscatter amplitude
datout(count,6)=mean(advdatg(:,6));     % burst averaged z beam backscatter amplitude
datout(count,7)=mean(advdatg(:,7));     % burst averaged x beam correlation coefficient
datout(count,8)=mean(advdatg(:,8));    % burst averaged y beam correlation coefficient
datout(count,9)=mean(advdatg(:,9));    % burst averaged z beam correlation coefficient
datout(count,10)=mean(advdatg(:,10));   % burst averaged heading
datout(count,11)=mean(advdatg(:,11));   % burst averaged pitch
datout(count,12)=mean(advdatg(:,12));   % burst averaged roll
datout(count,13)=mean(advdatg(:,13));       % burst averaged Temp (C) 
datout(count,14)=mean(advdatg(:,14));       % burst averaged pressure (counts)--uncalibrated
datout(count,15)=compass_head;
datout(count,16)=cv(2,2);             % u'u'
datout(count,17)=cv(2,3);             % u'v'
datout(count,18)=cv(2,4);             % u'w'
datout(count,19)=cv(3,3);             % v'v'
datout(count,20)=cv(3,4);             % v'w'
datout(count,21)=cv(4,4);             % w'w'
datout(count,22)=mean(Cb);            % burst averaged TSS concentration
datout(count,23)=cv(1,2);             % C'u'
datout(count,24)=cv(1,3);             % C'v'
datout(count,25)=cv(1,4);             % C'w'
datout(count,26:27)=b;
datout(count,28)=nsmp;        % number of samples remaining after quality control and used for averaging 
datout(count,29)=std(ur2);   % burst averaged U velocity standard deviation
datout(count,30)=std(vr);   % burst averaged V velocity standard deviation
datout(count,31)=std(wr);   % burst averaged W velocity standard deviation

datout(count,32)=Wspump;    % Ws  directly observed settling velocity using backCONC from TSS pump samples (mm/s)
datout(count,33)=Re_stress; %  u'w' directly observed turbulent Reynolds stress in (cm/s)^2
datout(count,34)=tau_cov;   % Directly observed stress in Pascals COV methodT

datout(count,37)=backConc2;                               %minimum concentration inferred from ADV backscatter

datout(count,38)=Ws_ALL;                                   %get settling velocity with background concentration included. Newly added 8/11/14.
datout(count,39)=mean_velocity;                           %save mean rotated U (KAF 2014)


datout(count,40)=Wsobs;         % Ws  directly observed settling velocity; back conc is the min CB observed over each depth per burst


%get TKE and TKE stress 
%TKE=1/2(u'2+v'2+w'2);
%TKE_tau=rho*0.2.*TKE; %off of  Biron et al., 2004 and Pope et al., 2006
%0.2 is the constant found by the above papers
datout(count,35)=0.5.*(datout(count,16)+datout(count,19)+datout(count,21)); %TKE 
datout(count,36)=P*0.2.*datout(count,35); %TKE derived stress
datout(count,39)=mean_velocity; %save mean rotated U (KAF 2014)

else
datout(count,:)=nan;


end
%datout(c,:) =[ U  V  W xamp yamp zamp  xcor ycor zcor headAvg pitchAvg rollAvg T dep  compass_head upr2 uvpr uwpr vpr2 vwpr wpr2 C Cupr Cvpr Cwpr  b     nsmp Ustd Vstd Wstd ws uw2 stress_obs];
 %              1  2  3   4  5    6     7    8    9     10     11       12      13  14   15          16   17    18   19   20   21  22 23  24   25   26 27 28    29   30   31  32  33  34
 end
 
 %eval(['save ' filename 'AVG_adjust_for_backConc.txt datout -ascii'])
 eval(['save ' filename '_AVG.mat datout bscatter_stats -mat']) %save as matlab workspace
 
end 